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Impact of dimensionless numbers on the efficiency of 
MRI-induced turbulent transport. 
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ABSTRACT 

The magneto-rotational instability is presently the most promising source of turbulent 
transport in accretion disks. However, some important issues still need to be addressed 
to quantify the role of MRI in disks; in particular no systematic investigation of the role 
of the physical dimensionless parameters of the problem on the dimensionless transport 
has been undertaken yet. For completeness, we first generalize existing results on the 
marginal stability limit in presence of both viscous and resistive dissipation, exhibit 
simple scalings for all relevant limits, and give them a physical interpretation. We 
then reexamine the question of transport efficiency through numerical simulations in 
the simplest setting of a local, unstratified shearing box, with the help of a pseudo 
spectral incompressible 3D code; viscosity and resistivity are explicitly accounted for. 
We focus on the effect of the dimensionless magnetic field strength, the Reynolds 
number, and the magnetic Prandtl number. First, we complete existing investigations 
on the field strength dependence by showing that the transport in high magnetic 
pressure disks close to marginal stability is highly time-dependent and surprisingly 
efficient. Second, we bring to light a significant dependence of the global transport 
on the magnetic Prandtl number, with a oc Pm s for the explored range: 0.12 < 
Pm < 8 and 200 < Re < 6400 (5 being in the range 0.25 to 0.5). We show that 
the dimensionless transport is not correlated to the dimensionless linear growth rate, 
contrarily to a largely held expectation. For large enough Reynolds numbers, one 
would expect the reported Prandtl number scaling of the transport should saturate, 
but such a saturation is out of reach of the present generation of supercomputers. 
Understanding this saturation process is nevertheless quite critical to accretion disk 
transport theory, as the magnetic Prandtl number Pm is expected to vary by many 
orders of magnitude between the various classes of disks, from Pm <C 1 in YSO disks 
to Pm > or 3> 1 in AGN disks. More generally, these results stress the need to control 
dissipation processes in astrophysical simulations. 



INTRODUCTION 



Angular momentum transport has always been a cen- 
tral issue in accretion d isk theory. The first a model 
|Shakura fc Sunvaevll 19731 ) already assumed the presence of 
strong turbulent motions, modelled through an effective vis- 
cosity, orders of magnitude larger than the expected disk 
molecular viscosity. Since then, the physical origin of this 
turbulence has been highly debated. As purely hydrody- 
namic non stratified Keplerian flows are known to be lin- 
early stable to small perturbations, a finite amplitude in- 
stability was first envisioned to trigger turbulence. This 
quest i on was studied both experimentally (|Richard fc Zahnl 
Il999l; lRichar3 l200ll ; iRichard et al]|200ll ) and numerically 
(jBalbus et al.l ll996: Ha wlev et al.1 '1999). leading to a long 
controversy. More recent numerical and experimental inves- 
tigations of this problem strongly support the idea that the 
transport due to this mechanism, if present, would be far to 
inefficient to account for the short disk evolution time-scales 



impo s ed by astroph ysical observations (|Lesur fc Longarettil 
120051 : I Ji et alj [2PO6) . Linear instabilities of hydrodynamic 
origin have also been envisioned as a source of tur- 
bulence, relating in particu l ar to the flow stratification 
felahr fc Bodenheimerl 120031; lUrpinl 120031 ; iDubrulle etafl 
120051 ; IShalvbkov fc Ruedigerl 120051 ). but these are ei- 
ther not present o r too i nefficient dJohnson fc Gammiel 
20061; lArlt fc Urpinl |2004| ; lLongaretti fc Lesurl 120071 : see 



Lesur fc LongaretJ 2 ^Tln d references therein fo~ecent 



review of this issue). 

The potential role of MHD instabilities in accre- 
tion disks was pointe d out in a seminal paper by 
iBalbus fc Hawlevl i|l99ll ), devoted to an analysis of what 
is now known as the magneto-rotational instability (MRI). 
This instability has extensively been studied since then , 
mainly with th e help of l ocal l[ Hawle y et al.l 1 19951 ; 
IStone et all Il996j ) and global (|Hawlevl |2000| ) 3D numeri- 
cal simulations. Although a more recent set of numerical 
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simu lations did focus on MRI energetics l|Gardiner fc Stone] 
l2005h . the dissipation of turbulent fields in these simulations 
is not controlled, as no physical term was introduced to ac- 
coun t for physical viscosity a nd resistivity. Note however 
that iBrandenburg et al.l (| 19951 ) have introduced such dissi- 
pation in their simulations, but kept it as small as possible, 
and in any case did not try to investigate their effect in a 
systematic way. Resistivity effects alone have also been in- 
troduced by iFleming et all (|2000t) . but viscous effects were 
still neglected. This raises questions about the exact role of 
numerical dissipation in all these simulations, especially at 
the light of our rec ent investigation of subcrit ical turbulence 
in accretion disks l|Lesur fc Longarettill2005l ). which clearly 
showed that a careful control of dissipation and resolution 
— and more generally of the dimensionless parameters of 
the problem — is required to properly quantify turbulent 
transport. 

This issue is addressed here in the context of MRI- 
driven turbulence, using a 3D spectral Fourier code, which 
allows a precise monitoring of viscous, resistive and numer- 
ical dissipation. First, we r ecall the MHD equ ations in the 
shearing sheet framework (lHawlev et al.|[T995l ). along with 
the relevant dimensionless parameters of the problem, and 
summarize what is known about their effect on MRl-induced 
turbulent transport. Next, we investigate the linear stability 
of the MRI, which to the best of our knowledge has not been 
fully characterized when both viscosity and resistivity are 
accounted for in the dispersion relation. Then, we present 
new results on the behavior of turbulent transport in dimen- 
sionless parameters regime that have not been investigated 
in the past: first, very close to the threshold of instability 
(in terms of relative magnetic field strength) , and then with 
respect to the magnetic Prandtl number, which has been ig- 
nored in all previous investigations. The dependence of tur- 
bulence transport on the magnetic Prandtl number is the 
most significant finding of this investigation. This depen- 
dence may turn out to be critical, as the magnetic Prandtl 
number varies by many order of magnitudes in astrophysi- 
cal disks. The astrophysical implications of our findings are 
further discussed in our concluding section, along with po- 
tential caveats relating to numerical limitations that may 
influence our results. 



2 SHEARING BOX CHARACTERIZATION 
AND SUMMARY OF EARLIER RESULTS 

The MRI has alr eady been ex tensively studied in the lit- 
erature (see, e.g., iBalbus! 12003! and references therein for a 
review of the subject). Our objective is to extend previous 
work through a systematic exploration of the dependence of 
the MRI-induced transport on the physical quantities char- 
acterizing the pro blem. For simplicity , we work in a shearing 
sheet setting (see lHawlev et alTll995l for a description of the 
shearing box equations, numerical boundary conditions, and 
conserved quantities); vertical stratification is ignored, but 
both viscous and resistive microphysical (molecular) dissipa- 
tion are included. This differs from previous investigations, 
where this is always ignored. Our previous experience with 
subcritical hydrodynamic transport has shown us that the 
inclusion of explicit dissipation is required to precisely char- 
acterize transport properties and ro sort out converged sim- 



ulatio ns from under-resolved ones (see iLesur fc Longarettil 
120051 for an extensive discussion and illustration of these 
points) . 

The problem is formulated in a cartesian frame centered 
at r = Ro, rotating with the disk at Q = Q(Ro) with radial 
dimension H <C Ro - In this work, H is the size of our simula- 
tion boxes, in the vertical and radial dimensions. This leads 
to the following set of equation, assuming <f> — *• x, r — » — y: 



d t u + u ■ Vu — 

dtB = 

V • 11 = 

V B = 



--VP+ — (V x B) x B, (1) 

P Hop 

— 2fl x u - 2Q.Sye v + uAu 

V x (u x B) + r]AB, (2) 

0, (3) 

0, (4) 



where the medium is defined by S — —rd r Q. For simplic- 
ity, incompressible motions are assumed. This is a priori 
justified by the fact that MRI-induced motions are usually 
subsonic, so that one expects at least in first approxima- 
tion that compressibility effects do not play a major role 
in the problem. This approximation allows us to remove the 
flow Mach number from the list of dimensionless parameters 
characterizing the problem, so that we can more effectively 
isolate and quantify the role of the various physical agents. 

The terms on the right-hand side member of Eq. ([TJ are 
the gas pressure, Lorentz force, Coriolis force, tidal force, 
and viscous dissipation, respectively. The steady-state so- 
lution to this equation is the local profile u = Sye x with 
S — 3/2Q, for Keplerian disks. In the remainder of this paper, 
we will use the deviation from the laminar profile w defined 
by w — u — Sye x . For simplicity, we assume that the steady- 
state magnetic field -Bo lies along the vertical axis. Note that 
this field is also the average field in the shearing sheet box, 
and is conserved during the evolution thanks to the shearing 
sheet boundary conditions l|Hawlev et al.|[l995l ). 

These equations are characterized by four dimensionless 
numbers, the first three relating to the Navier-Stokes equa- 
tion while the last one belongs to the induction equation: 

• The Reynolds number, Re = SH 2 jv, measuring the 
relative importance of nonlinear coupling through the ad- 
vection term, and viscous dissipation. 

• A proxy to the plasma beta parameter, defined here as 
3 = S 2 H 2 /Vl where V% = B 2 /fi p is the Alfven speed. The 
rationale of this definition follows from the vertical hydro- 
static equilibrium constraint c s ~ Q.H, which is expected 
to hold in thin disks, so that our definition of 3 is indeed 
related to the plasma parameter in an equivalent, vertically 
stratified disk. This parameter measures the relative weight 
of the Lorentz force and the advection term. 

• The rotation number (inverse Rossby number), defined 
as Rn — 20./ S, which measures the relative importance of 
the Coriolis force. 

• The magnetic Reynolds number, R m = SH 2 /rj, which 
measures the relative importance of resistive dissipation 
with respect to the ideal term in the induction equation. 

We consider only Keplerian disks in this investigation, 
so that the rotation number is held fixed to its Keplerian 
value Rn = —4/3. This leaves us with three independent 
dimensionless numbers: 8, Re, and Rm. 
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On the other hand, the (local in the disk) dimensionless 
transport coefficient, 



(v x Vy - B X By/(fXop)) 



(5) 



S 2 R 2 

being a dimensionless number, can only depend on the lo- 
cal dimensionless parameters characterizing the flow that 
we have just defineo{3 (the bracketing refers to appropri- 
ate box and/or time averages). Our task reduces to char- 
acterize this dimensionless transport, as a function of the 
three independent dimensionless numbers just defined. How- 
ever, for later convenience, we take them to be (5, Re and 
P m = v/r] = Rm/Re instead (the rationale of this latter 
choice will become apparent later on). 

All previous investigations ignore the dependence on 
the last two dimensionless numbers, who have not been in- 
cluded in the physi cal description up t o now. Within such 
an approximation, lHawlev et ail l| 19951) have characterized 
the dependence of a on j3. Their results imply that 



3/3 



-1/2 



(6) 



from their Eqs. (10), (15), (16) and (18). 

This implies in particular that a increases when the ini- 
tial (and box average) magnetic field B is increased. How- 
ever, for a large enough field, the smallest unstable wave- 
length (which increases along with B ) becomes larger than 
the box size, and the instability is quenched. On this ba- 
sis, one expects that the scaling Eq. ((6]l would break down 
close enough to the critical [3 stability limit. This question is 
somewhat investigated in the present work. However, most 
of our effort is devoted to characterizing the Re and Pm 
dependence of a. 



3 LINEAR STABILITY ANALYSIS 

In order to quantify the MRI induced turbulent transport, 
it is first necessary to define the parameter domain in which 
this instability operates. The linear stability of differentially 
rotating disks in presence of a magnetic field was first in- 
vestig ated in the astrophysical context bv lBalbus fc Hawlevl 
|l99ll ). Then, the in stability in the weakly ionized case 
has been considered dBlaes fc Balbusl Il994l ; IWardlei 1 19991 ; 
iBalbus fc Terquemll200ll), lea ding to the well known Dead 
Zone problem l|Gammie 1 19961 ). However, we are not aware 
of any reasonably complete and heuristically clarified inves- 
tigation of the stability limits of the fluid when both viscous 
and resistive dissipation are taken into account. Some dis- 
cussions of this point are available in the literature, mostly 
motivated by l i quid-metal experiments, in t he limit Pm <C 1 
jji et al.ll200ll ; iRiidiger fc Shalvbkovll2002l) . However, these 
papers exhibit no clear asymptotical limit that may be use- 
ful for astrophysical disks. Therefore, we provide such an 
analysis here, as a prelude to our nonlinear simulations. 

We will consider only axisymmetric perturbations, so 
that we can eliminate the azimuthal perturbation transport 
term. Note that this assumption does not seem to have a 



great influence on the stability limit, since 3D numerical 
simulations and linear analysis of axisymmetric modes ex- 
hibit nearly the same stability limit; this holds in particular 
in the simulations presented here. 

We linearize and Fourier transform the equations of 
motion by assuming v = t)oexp(j(tjt — k y y — k z z)j and 
b — bo exp (i(ut — k y y — fc z z)). This yields the following 
linearized equation set: 



(ico + vk 2 )vo = i kip — ik z ^° bo 

+ (2f2 — S)v y e x — 2Q.v x e y 
(iuj + rik 2 )b = -ik z B v + b y Se x , 
ik ■ v — 0, 
ik B = 0, 



(8) 
(9) 
(10) 



where ip is the perturbation in total pressure (P + B 2 / /j,o) / p. 
Introducing u> v = lu — ivk 2 and lu v = to — ir/k 2 , the Alfven 
speed V\ = B 2 /fj,op, the epicyclic frequency k 2 = 2fi(2f2 — 
S) and j 2 = k 2 /k 2 , one eventually gets the dispersion rela- 
tion: 



{iOuLOri — k 2 z V]i) I LU 2 LU 2 — 2iU l ,LU v k 2 VA ~ U) 2 ^ '' 



-k 2 z Vl{2Q,S-i 2 - k 2 z Vl^j = (11) 



which we now solve in various dissipation regimes. 



3.1 Pm = 1 behavior 

Let us first look at the Pm — 1 case, where the dispersion 
equation can be solved exactly by analytical means. The 
condition Sr(w) < expresses the existence of the instability, 
and implies that the MRI exists if and only if u 2 k 4 < —lu 2 . 
From this constraint and the dispersion relation Eq. {TTJ, 
one finds that: 



It may also depend on the simulation aspect ratio and resolu- 
tion, from a numerical point of view. 



2 + 16k 2 V 2 Q2j 2 k 2y2 K 2^2 

V < 2fc* fc4 2k 4 > 

is a necessary and sufficient criterion for instability. One 
can check that the highest v values obtain when 7=1 and 
k z — min(fc z ) = 2-n/H, which corresponds to the so-called 
channel flow solution in the z direction. From our definition 
of the Reynolds number as Re — SH 2 /u where H is the 
numerical box height or the typical disk height, and of the 
plasma parameter (3 — S 2 H 2 /V%, the stability limit Eq. {12} 
translates into a relation between these two parameters, rep- 
resented on Fig. [T] 

Note that the instability has two different limits, de- 
pending on the P parameter: 

• A high f3 regime, corresponding to a low magnetic pres- 
sure. In this regime, marginal stability occurs at a character- 
istic Reynolds number value Re c — 80. This behavior illus- 
trates that the growth time scale of the most unstable mode 
must be shorter than the dissipation time scale, defined by 
T d ~ k 2 jv. 

• A low f3 regime, which is nearly Reynolds independent. 
In this region, one can define a critical f3 (/3 C = 29.5) for 
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Figure 1. MRI linear stability limit for Pm = 1 



which the MRI is lost. This behavior can be explained by 
considering the unstable mode of shortest wavelength: as f3 
goes to smaller values, the smallest unstable wavelength in- 
creases (see Eg. I12[) . At some point it becomes larger than 
the scale height H (or box size in our case) and the insta- 
bility is lost. Since this phenomenon takes place at large 
scales, the Reynolds number plays little or no role. Note 
that this regime is not specific to our unstratified calcula- 
tion, since similar results are found for a stratified medium 
where marginal stabilit y us ually occurs for /3 e ^ 1 ( see e.g 
iBalbus fc Hawlevlll99ll and lGammie fc BalbusHl994[ ). This 
limit is reached when the last factor in Eq. cancels out, 
i.e., when 2Q.S = V%k 2 z (the usual dissipationless MRI sta- 
bility limit). 



3.2.1 Low viscosity limit: 

First consider the limit where vk 2 <C k. In this case, Eq. (|13|l 
reduces to r\k 2 = (2QS '/ 'k 2 ) 1 ' 2 !^^ (where rjk 2 <C k has 
been self-consistently used). Using the Lundquist number 
defined as Lu = iJm/J -1 ' 2 , this can be recast as 

^=(J^)-3.6. ("J 

Note that our defi nition of the Lundq uist number is not 
strictly identical to iTurner et all (|200ft ) but is widely used 
in the MHD communitjjf] In this regime, the lo 2 k 2 term bal- 
ances the 2Q.SV\k\ term in the dispersion rel ation Eq . (1111) . 
Eq. 11 141) corresponds to the limit found by iFleming et al.l 
(2000). It is also related to the origin o f the "dead zone" in 
accretion disks (see e.g iGammi el ll996l ). This marginal sta- 
bility limit is relevant to disks with low Prandtl numbers 
(Pm 1) and high Reynolds numbers (Re 3> 1), such as 
YSO disks. 

Also, for negligible resistivity, growth rates in this 
regime are given by 

This result is valid for Va^z < k due to our expansion 
scheme; it also gives the correct order of magnitude of maxi- 
mum growth rates when Vkfcz shown by the standard 
dissipationless MRI analysis. 

3.2.2 High viscosity limit: 

Conversely, consider the large viscosity limit, where vk 2 
k. The corresponding relations in this limit are 



3.2 Pm 7^ 1 behavior 

The dispersion relation can no longer be solved exactly in 
this case, but an approximate solution can be found in 
the low magnetic field limit (Va — > 0, or more precisely 
Vkfcz <C ft), where marginal stability follows from a bal- 
ance between the destabilizing term, and the dissipation 
ones. The "opposite" (high f3) marginal stability limit, where 
destabilization is balanced by the usual Alfvenic magnetic 
tension, is briefly addressed at the end of this section. 

In the limit of vanishing magnetic field, the dispersion 
relation has two relevant roots ui 2 — and u> 2 = k 2 . In 
what follows, we refer to these roots as the Alfvenic and the 
inertial branch, respectively. Looking for the first order cor- 
rection in V%k 2 to the Alfvenic branch yields the following 
result, which describes the MRI modes: 

U = ivk ° ± \ k}{ V - n J)* + K * ) VAk - (13) 

Note that viscosity and resistivity do not play a sym- 
metric role in this expression. Two interesting limits with 
respect to the magnitude of the viscosity prove useful to 
characterize marginal stability. As before, we maximize in- 
stability by assuming 7 = 1 and k z = 2n/H . 



nl/2 

ReRm = ^-(2tt) 3 /3 1/2 ~ 215/3 1/2 . (16) 

and 

_i 1 /2(IS\ Tr , 

t ^fer^- (i7) 

In this regime, the u) 2 u> 2 term balances the 2Q,SV\k\ term 
in the dispersion relation Eq. (|11[) . The growth rates rele- 
vant here are much smaller than in the small viscosity limit, 
Eq. (|15|) . In fact, Eq. (|13[) indicates that this is the case as 
soon as vk 2 < k, or equivalently, for the largest mode, when 

Re > 3(2tt) 2 /2 ~ 60. (18) 

This limit divides the low and high viscosity regime. 

The marginal stability limit Eq. (|16[) obtains for large 
Prandtl and small Reynolds numbers. In the large Prandtl 
(Pm S> 1) and large Reynolds number limit (Re 3> 1) ex- 
pected in AGN disks, the growth rates of Eq. (|15[) , or more 
generally of dissipationless MRI, are recovered. As before, 

2 The difference lies in the fact that our calculation is made in 
the limit of high /3, leading to a linear growth rate controled by 
Va instead of Q 
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these growth rates are expected to be valid (in order of mag- 
nitude) for Vkfez < « due to our expansion scheme. 

Note finally that a similar analysis can be performed 
for the inertial modes, but is not very informative; as they 
appear to be always stable. 



3.2.3 High f3 limit: 

Although we did not investigate this case in much detail, 
it is apparent from Eq. that when 2QS = V\k\ [can- 
cellation of the last term in Eq. pip ]. u v = is one of 
the solutions to the dispersion relation. At the light of our 
preceding analyzes, and because this equality embodies the 
MRI stability limit in the ideal case, as recalled above, it is 
apparent that this relation is the relevant limit in a small 
dissipation context as well, generalizing the result found for 
Pm = 1. 



3.2.4 Heuristic explanation: 

To explain the behavior brought to light in Eqs. (|14[) and 
(|16p . it is useful to recall the physical origin of th e insta- 
bility, as discussed, e.g., in lBalbus fc Hawlevl (|2003f ). in the 
dissipation-free limit; the process is sketched on Fig. [2] for 
convenience. Assume for definiteness that one starts by dis- 
torting the equilibrium velocity field in the radial direc- 
tion with a sinusoidal perturbation in the vertical direction: 
v y — v yo exp ( — ikz) . The magnetic field being frozen in 
the fluid will also develop a radial component [first term 
in the right-hand side member of the linearized induction 
equation, Eq. ||B)]: the shear will then transform this radial 
field in an azimuthal one [second term in the right-hand 
side member of the linearized induction equation, Eq. JHJ]. 
The resulting tension force produces a momentum trans- 
fer between fluid particles that have been moved according 
to the imposed velocity perturbation [second term in the 
right-hand side member of the linearized motion equation, 
Eq. @], This force is destabilizing if the angular velocity 
decreases with radius: indeed in this case, the inner parti- 
cle, moving faster than the outer one, will transfer orbital 
momentum to the outer one, thereby reinforcing its inward 
motion, an effect mediated by the Coriolis force when seen 
in the rotating frame. In this description, marginal stabil- 
ity follows when the driving mechanism is balanced by the 
usual tension restoring force (the piece not connected to the 
generation of magnetic field from the mean shear). 

What does dissipation change to this picture ? For defi- 
niteness, let us focus on marginal stability and let us consider 
only resistive dissipation for the time being ( "low" viscosity 
limit). In this limit, the magnitude of the velocity and mag- 
netic fields in the various steps of the instability mechanism 
described above are controlled by dissipation processes so 
that one may again go through the preceding process step 
by step, assuming equilibrium at each step. The magnitude 
of the radial magnetic field in this context results from the 
balance between the motion driving and field dissipation: 



ikBoVy = nk b y 



(19) 



while the shearing generation of the azimuthal field from the 
radial one is also balanced by resistive dissipation: 



dr 2 n/dr > 
dn/dr < 
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Figure 2. Sketch of the MRI mechanism (see text). 



Sb v 



rjk 2 b x . 



(20) 



Both relations follow from the induction equation in the 
marginal stability limit, except for the term dropped in 
Eq. (|20[) , which leads to the usual magnetic tension sta- 
bilization and is of no interest in the limit considered here. 
The azimuthal force balance then requires that 



(20 - S)v y 



I O x , 

fJ-oPo 



(21) 



* — 2£lSVj±k , once the two preceding constraints 
are taken into account (inclusion of u> in this line of argument 
does not change the result). As noted earlier, this relation 
directly leads to Eq. 

If one assumes instead that viscous dissipation exceeds 
the Coriolis force in magnitude, then the magnetic tension 
due to the generation of azimuthal field from the radial 
one by the shear should be balanced by viscous dissipation 
instead of the Coriolis force in the two horizontal compo- 
nents of the momentum equation, leading alternatively to 
= 20SVjfc 2 , i.e. to Eq. (JTSJ. 
This also relates to the structure of MRI modes. In the 
limit of a very small magnetic tension restoring force, the 
Alfvenic branch is made of b x dominated modes. The other 
components of the magnetic field and the velocity field are 
of the order of V/tfc compared to b x . Therefore, the growth 
rate is at first controlled by the dissipation rate of b x , which 
is related to the resistivity [first term of the right hand side 
of Eq. Q3] • The interaction of the other fields, which leads 
to the MRI, is controlled by a term symmetric in v and r\ 
[second term of Eq. I13[l ] . 
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Figure 3. MRI linear stability limit in the Pm ^ 1 case for 
P = 10 4 . 



Figure 4. MRI growth rate (arbitrary unit) as a function of 
viscous and resistive dissipation for j3 = 10 4 



3.2.5 Generic behavior: 

A more complete view of the stability limits and growth 
rates implied by Eq. may be obtained from exact nu- 
merical solutions for Pm 7^ 1. Expressing this dispersion 
relation in terms of u> leads to the condition: 



u -2ik 2 Lj 3 (n + v) 



• fcV 



+ 4rjv) + b 



+uj[2ik 6 (r)V 2 + vrj 2 ) + aik 2 {v + ri) + 2ibi]k : 



+v 2 r) 2 k 8 + avnk 4 + bn 2 k* 



0, 



with 

a — 
b = 
c = 



2k\v\ 
2 2 

K 7 



k 2 z VX(2QYS- 



klvl) 



(22) 

(23) 
(24) 
(25) 



To characterize the stability limits as a function of 
the Reynolds and the Magnetic Reynolds number {Rm = 
SH 2 /rj), one needs to choose j3, 7 and k z . As in the Pm = 1 
case, we take k z — 2-k/H and 7 = 1 (which are again ex- 
pected to maximize the dissipation limits), and solve the 
relation (I22|l for j3 = 10 4 . The resulting stability limits 
are shown on Fig. [3] and the corresponding growth rates 
on Fig. [4] (arbitrary units). These results match closely 
the analytical limits just discussed: a high Re threshold 
found for Rm ~ 371, and a low Re threshold found for 
RmRe ~ 2.3 x 10 4 , both in agreement with Eqs. (|14[) and 
(|16p . respectively. Moreover, significantly lowered growth 
rates are observed when Re 60 to 80, as predicted by 
Eqs. (|17fl and (|18[) . A similar behavior follows at much 
smaller j3. For example, the observed scalings are identi- 
cal, and the preceding asymptotic expressions valid within a 
factor of 2, down to j3 values of the order of twice the critical 
f3 limit. 

These results indicate that most of the stability limit 
behavior is captured by the approximate relations Eqs. (|14p 
and (|16|) (as well as by the large field (5 limit, where relevant) , 
whose physical origin has been discussed above. 



3.3 Numerics 

3. 3. 1 Equations 

Our objective is to simulate the system of Eqs. (JXJ) and ([2]), 
with the incompressiblity condition Eq. (|3}, to characterize 
the dependence of turbulent transport on the main dimen- 
sionless numbers introduced above (/3, Re and Pm). We 
focus on incompressible motions; indeed, the values of a 
found in previous investigations makes us a priori expect 
that compressibility effects will be small. In any case, this 
allows us to more effectively distinguish the effects of the 
various physical mechanisms at work in this problem. 

First, we simplify the problem from a numerical point 
of view by distinguishing the mean laminar shear u = Sye x 
from the deviation from this mean w. The resulting equa- 
tions read: 

<-r <-, , B x (V x B) 

dtw + w-vw = — byo x w — vip H 

+ (2fi — S)w y e x — 2Qw x e y + vAw 
-Syd x B + B ■ Vw + B y Se x + nAB 



d t B + w ■ VB 
V • w 



= 



V • B = 

This system is numerically solved using a full 3D spec- 
tral cod e, using the classica l shearing sheet boundary con- 
ditions l|Hawlev et al.l 1 19951 ) . This code is now briefly de- 
scribed. 



3.3.2 Numerical code 

The code used for thes e simulations is an MHD extension 
of the HD code used in lLesur fc Longarettil (|2005l ). and ex- 
tensively described there. This code is a full 3D spectral 
(Fourier) code, based on FFTW libraries, parallelized using 
the MPI protocol. This kind of code has many advantages 
for the simulation of incompressible turbulence, such as: 

• The incompressibility and solenoid conditions are easily 
implemented at machine precision, using a projector func- 
tion in Fourier space. 
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• The energy budget is much easier to control, leading 
to a precise quantification of the energy losses by numerical 
dissipation. 

• Spatial derivatives are very accurate down to the grid 
scale (equivalent to an infinite order finite difference scheme 
down to the grid scale). 

The algorithm used is a classical pseudo spectral 
method which may be described as follows. All the deriva- 
tives are computed in Fourier space. However the nonlinear 
term require special treatment : in Fourier space, a real space 
product is a convolution, for which the computational time 
evolves as n 2 , where n is the number of grid cells. The com- 
putation time is minimized if one goes back to real space, 
compute the needed product and then transforms the re- 
sult to Fourier space. This procedure (pseudo spectral pro- 
cedure) is more efficient than a direct convolution product 
since the FFT computation time scales as n log n. However, 
the finite resolution used in this procedure generates a nu- 
merical artifact commonly known as the "aliasing" effect 
(apparition of non physical waves near the Nyquist Fre- 
quency). This effect may be handled through a dealiazing 
procedure, in which the nonlinear terms are computed with 
a resolution 3/2 highe r than the eff ective resolution used in 
the source terms (e.g.. |Pevretll2002h . 

Comparing our spectral code with a ZEUS-type finite 
difference code l|Stone fe Normanl [l992 j. similar results are 
obtained with a finite difference resolution two to three times 
higher than the spectral resolution. However, FFTs calcula- 
tions are more computationally expensive than finite differ- 
ences, leading to a final computational time equivalent for 
both kind of code with the same "effective" resolution. 

All the simulations presented in this paper were per- 
formed with an xyz resolution of 128 x 64 x 64 with an 
aspect ratio of 4 X 1 X 1, a; being the azimuthal direction, y 
the radial direction and z vertical direction. One may change 
the physical viscosity and resistivity as well as the magnetic 
field intensity (/3). The mean magnetic field (conserved in 
the simulations due to the adopted boundary conditions) is 
aligned in the z direction. White noise initial perturbations 
with respect to the laminar flow are introduced as initial 
conditions on all variables. With (3 = 100, Pm = 1 and 
Re = 3200 one typically generates flow snapshots as shown 
on Fig. [5] after relaxation of transients; this flow is quite 
characteristic of a fully developed 3D turbulent field 0. 

3.4 MRI behavior near the instability threshold 

The MRI is a weak magnetic field instability, which should 
be quenched for /3 < 1 in astrophysical disks. Since the 
MRI is assumed to be the source of momentum transport 
in disks, and as at least some disks are expected to be close 
to equip artition if they are to support magnetically driven 
ejection |Ferreiral Il997f ). on may wonder if this instability 
is efficient enough in the vicinity of the strong magnetic 
field stability threshold. We investigate this question in an 
unstratified context here (the absence of stratification sig- 
nificantly raises the f3 stability threshold). 

3 Movies of some of the simulations presented in this paper may 
be found on the web at 

http://www-laog.obs.ujf-grenoble.fr/public/glesur/index.htm 




Figure 5. w y plot (radial velocity) for /3 = 100, Re=3200, 



We present two simulations. In the first one, (3 — 100 
and Re — 3200 (run 1); this run reproduces typical results 
from the literature. The second simulation is performed close 
the /3 threshold, i.e. for (3 = 30 and Re = 3200. The time de- 
velopment of some important quantities is depicted on figs 
(0 and for these two runs. One immediately notices a 
marked difference between these two simulati ons. On run 1, 
we fin d a classical MRI behavior, as studied bv lHawlev et al.l 
(1995), characterized by a ~ 10 _1 and random fluctuations 
of all statistical quantities. However, run 2 exhibits strong 
exponential growth ("bursts") for about 100 shear times 
(~10 orbits), and a sudden drop of fluctuation amplitudes. 
This behavior is explained as follows: for such low (3 only 
the largest wavelength mode is unstable (and not smaller 
scales); the mode amplitude increases for many shear times, 
as this mode is an exact nonli near solution to the in com- 
pressible equations of motions (jGoodman fc Xulll994l l. We 
therefore observe t he gr owth of the channel flow as seen by 
lHawlev fc Balbusj l|l992ft . However, as this channel solution 
reaches sufficiently large amplitude, secondary instabilities 
such as the Kelvin-Helmoltz instability quickly set in and 
destroy these ordered motions, and a new cycle starts (see 
iGoodman fc Xull 19941 for a detailed description of these sec- 
ondary instabilities). 

Note that this kind of behavior and related explana- 
tion does also apply to the low Reynolds threshold, since 
there the smallest scales are viscously damped and only the 
largest ones remain unstable. We did obser ve this behav- 
ior clo se the low Reynolds threshold, as did iFleming et al.l 
(2000) but in an indirect way (see Figs. 2 and 4 of their 
paper ), and one can conclude that these bursts are charac- 
teristic of a marginally unstable MRI. Such bursts may be 
astrophysically relevant. Indeed, one may wonder ab out the 
MRI behavior close to the dead zone (|Gammielll996h . where 
presumably the magnetic Reynolds number is small, and the 
instability quenched. If these bursts exist in real disks, they 
may quickly destroy this dead zone under the effects of the 
strong turbulent motions observed in our simulations. 

Let us have a closer look on these bursts with the help 
of correlation lengths defined as 



U 



J dyi J f(xi)f(xi - yi)da 



(26) 



/ f 2 (xi) dxi 

where i — 1, 2, 3 is the direction of correlation and / refers ei- 




ther to the velocity or magnetic field. Note that with this def- 
inition, the correlation length vanishes in the z direction for 
a pure sinusoidal signal, and equals 1 in the y direction for a 
channel flow, as a consequence of the shearing sheet bound- 
ary conditions. Therefore, these correlation lengths provide 
us with a convenient tool to trace the presence of the channel 
flow solution in our simulations. We show on fig.[8]and fig. [9] 
the evolution of the correlation length in the y and z direc- 
tion for the w x field (a similar behavior is obtained with the 
other field components). The behavior of correlation lengths 
closely follows what can be seen by monitoring the energy 
in the deviations from the laminar flow (fig. 01, and indicate 
the presence of two main regimes in this simulation. The 
first regime corresponds to an exponential growth ("burst") 
of the channel flow, for which L y is found to be equal to 
the box size and L z = (a careful examination shows that 
L z is exponentially decaying down to 10 -10 ), indicating the 
presence of a purely sinusoidal mode in the z direction in 
the burst stage. The second regime is a more classical state 
for 3D turbulent motion, with L y ~ 0.5 and L z ~ 0.4. Note 
that L y grows on very short time-scales, leading eventually 
to a new burst stage. 

These correlation lengths disclose numerical artifacts in 
the first regime. In a real disk, one would expect a loss of cor- 
relation in the radial direction on a scale of the order of a few 
scale heights: indeed, the typical frequency involved in these 
phenomena is of the order of the Keplerian frequency and a 
signal can't propagate faster than the sound speed, leading 
to a maximum correlation length of a few scale heights. 

Similarly, the vanishingly small vertical correlation 
length for the channel flow solution is also an artifact of the 
adopted boundary conditions. A more realistic result would 
follow if one were to take into account the vertical strati- 
fication and set the boundary conditions far enough from 
the disk midplane. More generally, our results are probably 
not directly applicable to a real disk, but they shed some 
light on what the generic behavior of the MRI would look 
like near various stability thresholds, even though different 
aspect ratio and boundary conditions should be investigated 
before firm conclusions can be drawn. 

Finally, the behavior exemplified in our simulations sug- 
gests that assuming a constant would poorly represent the 
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Figure 8. Correlation length of w x in the y direction as a function 
of time, P = 30. 

transport behavior close enough to the marginal stability 
limit. Time-dependent transport models are needed in such 
a context. Real disks may not operate close to the strong 
field limit unless some (unknown) back-reaction loop is at 
work, or unless (more realistically) the magnetic field varies 
in a systematic way with radius throughout the disk; con- 
sequently, the bursting behavior observed here may imply 
a similar ejection variability in the relevant regions of jet- 
driving disks. Note however, that our "mean" equivalent a 
is rather large (a ~ 5), leading to question the role of the 
ignored fluid compressibility in these cases; it is quite pos- 
sible that couplings to compressible modes may effectively 
limit the magnitude of the bursts. 

3.5 Magnetic Prandtl dependence of transport 
coefficients 

All previously published simulations were performed with- 
out numerical control of the dissipation scales and dissipa- 
tion processes. However, as pointed out earlier, such a con- 
trol is required to ascertain convergence. In this section, the 
role of the Reynolds and Prandtl numbers defined in section 
[3] is examined. In particular, the Prandtl number allows us 
to change the ratio of the viscous and resistive dissipation 
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Figure 9. Correlation length of w x in the z direction as a function 
of time, /3 = 30. 



0.75 
0.7 

0.65 
0.6 

0.55 
0.5 













-0- Re=200 




—1— Re=400 




-•- Re=800 




-K-Re=1600 




-6- Re=3200 




Re=6400 



0.12 0.25 



1 

Pm 



Figure 11. Linear growth rate of the largest mode for various 
(Re, Pm) at /3 = 100 
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Figure 10. Prandtl effect for (3 = 100 



scales. Unfortunately, deviations from Pm = 1 are quite de- 
manding numerically, since one wants to resolve both the 
velocity and magnetic dissipations scales. We present on 
Fig. \W\ the result of such simulations: we plot the mean 
transport coefficient (a) as a function of the Prandtl num- 
ber, for various Reynolds numbers (the Reynolds number 
quantifies the viscous dissipation scale). Statistical averages 
are computed over 500 shear times, and start after the first 
100 shear times to avoid pollution by relaxation of the initial 
transient dynamics. From these plots, one finds a significant 
correlation between the Prandtl number and the transport 
coefficient, leading to 



q oc Pm 



for 



0.12 < Pm < 8 
200 < Re < 6400 



(27) 



with S in the range 0.25 — 0.5. Note that this results show 
that the transport coefficient depends on Re and Rm via 
Pm, at least in the Pm range concidered in this paper. This 
may be seen on Fig. [TO] as a small vertical dispersion (vari- 
ation of both Re and Rm at constant Pm) compared to 
the effect of a single Pm change. Although this section is 
the briefest of the paper, this result constitutes the most 
important finding of this investigation (and the most com- 
putationally intensive one!). 



Note that the numerical results obtained at very 
high Reynolds number and high Prandtl number are only 
marginally resolved, mainly because of a very short mag- 
netic dissipation scale. This remark may explain that the 
two points at Pm = 8 lie somewhat below the mean of 
the other results. Our preliminary tests at higher resolution 
seem to show that a higher transport obtains at higher reso- 
lution at Pm = 8 and Re — 6400, which confirms a limit due 
to resolution in these high Pm runs. This behavior is easily 
understood, since the finite numerical resolution enforces a 
numerical dissipation scale (roughly equal to the grid scale) , 
which is obviously the same for the magnetic and velocity 
fields. Therefore, at high Pm, the effective magnetic dis- 
sipation scale is forced to be larger than the expected one, 
leading to an altered spectral distribution and a smaller "nu- 
merical Prandtl". 

One may wonder if this Prandtl dependence may be 
correlated to the linear growth rates discussed before. To 
this effect, we plot the linear growth rate of the largest 
mode for the different simulations used for this study on 
Fig. 1111 Similar results follow when replacing the growth 
rate of the largest mode by the maximum growth rate. Al- 
though the idea of a transport efficiency controlled by the 
linear growth rate is widely spread in the Astrophysical com- 
munity, this plot shows us that, at least for this example, 
the linear growth rate doesn't explain the transport behavior 
observed on Fig. 1101 Moreover, it appears that, as one may 
suspect from Eq. (|13[1 , the growth rate is not controlled only 
by Pm, but also by some complex combination of Re and 
Pm. lUmurhan et all (|2007h tried to get this kind of a — Pm 
correlation analytically, using a weakly non linear analysis 
of the channel flow. This study leads to a stronger a — Pm 
correlation with S — 1 in the limit Pm <C 1, which appears 
to be quite different from our full 3D numerical results. Note 
however that their analysis belongs to very different bound- 
ary conditions (rigid instead of shearing sheet) and there 
results are therefore not of direct relevance to our numerical 
investigation; nevertheless both studies point out the role 
of the Prandtl number. In any case, one needs to find some 
full nonlinear theory to explain the transport dependance 
on Pm. 

The observed correlation indicates the existence of a 
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back-reaction of the small magnetic field scales on the large 
ones (at least for the range of Reynolds and Prandtl num- 
bers explored here). Note that this effect is expected to 
saturate at some point, since in the limit Pm — > with 
Re — > oo and Rm kept constant, Eq. (I27f) predicts a van- 
ishing transport in spite of the existence of the linear in- 
stability. This issue is further discussed in the conclusion, 
the Reynolds number limitation of our investigation being 
the most serious here. In any case, the exact implications 
of these findings remain to be understood, but may poten- 
tially be quite important since the Prandtl number varies 
by many orders of magnitude in astrophysica l objects. For 
example, iBrandenburg fc Subramanianl (|2005l ) suggest that 
values as small as Pm ~ 10 -8 might be found in young stel- 
lar objects, while Pm ~ 10 4 would be more typical of AGN 
disks. These estimates are highly uncertain; even a substan- 
tially narrower range is of course out of reach of present day 
computers. 

Finally, this kind of back-reaction points out the po- 
tential role of small scale physics (dissipation scales) on the 
properties of turbulence at the largest available scales (disk 
height scale). This argues for a careful treatment of the role 
of dissipation and reconnection processes on the turbulence 
transport characterization. 



4 DISCUSSION 

In this paper, we have investigated the role of local di- 
mensionless numbers on the efficiency of the dimension- 
less turbulent transport. To this effect, we have first gen- 
eralized previously published linear stability limits, to ac- 
count for the presence of both viscous and resistive dis- 
sipation. Namely, we have confirmed in all cases that the 
large field marginal stability limit is characterized by a 
constant plasma /3 parameter, of order 30 in the shearing 
sheet unstratified context (but more likely of order unity 
in real, stratified disks). When marginal stability follows 
from dissipation and not magnetic tension stabilization, wc 
have found that the marginal stability limit is captured by 
two asymptotic regime: a large Reynolds (Re), small mag- 
netic Reynolds one (Rm), with a marginal stability limit 
Rm ~ /3 1 / 2 , and a small Reynolds, large magnetic Reynolds 
number one, where ReRm ~ lO 2 /? 1 / 2 . A phenomenological 
explanation has been provided for this behavior. 

In the previous section, we have investigated the be- 
havior of the MRI near the low j3 instability threshold; in 
our simulations, (3 = 30, a value representative of the large 
field threshold in our simulation box. In vertically strati- 
fied disks, this threshold obtains for much smaller values, 
typically j3 ~ 1 (Gammie & Balbus 1994). We found, some- 
what unexpectedly, that turbulent transport is significantly 
enhanced through burst events, even surprisingly close to 
the marginal stability threshold. As pointed out earlier, this 
behavior is physical and not numerical. The use of peri- 
odic boundary conditions (vertical) or semi periodic (radial) 
boundary conditions may enhance the role of the channel 
flow solution which is responsibly for this behavior, and a 
real disk channel flow may break up sooner than observed in 
our local simulations, leading to smaller burst magnitudes. 
Moreover, ct > 1 leads to supersonic motions and compress- 
ible numerical simulations are needed to properly quantify 



the phenomenon, which may exhibit new secondary com- 
pressible instabilities in such a context. All these issues lead 
to the conclusion that low /3 MRI would produce weaker 
bursts and therefore smaller transport coefficient than ob- 
served in our simulation. However, there is no physical rea- 
son why the turbulence bursts would be suppressed, and 
we believe that these bursts may be a strong signature of 
regions of stratified disks where MRI-driven turbulence is 
driven close to the marginal stability threshold. 

The most important new result reported in this paper 
is a correlation between the transport efficiency, and the 
magnetic Prandtl number, leading to a higher transport co- 
efficient for larger Prandtl numbers. As in the case of the 
bursting behavior discussed above, the boundary conditions 
used in these simulations play some role in the result. How- 
ever, the possible biases are less obvious and tests with plane 
radial walls need to be performed to get a grasp on bound- 
ary condition effects. Moreover, one needs to check the cor- 
relation at higher resolutions, and if possible higher Prandtl 
numbers, using different kind of codes to get a better char- 
acterization and a physical description of the phenomena 
involved in this observation. 

More specifically, a puzzling fact points towards a po- 
tential bias due to the shearing sheet boundary conditions. 
In non-magnetized shear flows, transport in the subcritical 
regime, far enough from the marginal stability limit scales 
like 1/R a where Rq is the subcri tical transition Reynolds 
number ( Lesur fc Longaretti| [2005'l. Closer to the marginal 
stability limit, and in the supercritical regime (e.g., when 
the Rayleigh stability criterion is not satisfied), transport 
is enhanced with respect to this scaling, but one always has 
a < 1/ Rc where Re is the critical Reynolds number of linear 
instability. However, for MRI-driven turbulence, one has a > 
or^> 1/ Rc, as can be checked from our results. Close to the 
marginal stability limit, this enhanced efficiency is related 
to the existence of the channel flow solution, as discussed 
above. As each linear mode is a nonlinear solution to the in- 
compressible problem, one may ask whether this enhanced 
transport, which is observed also far from the marginal sta- 
bility limit, is not an artifact of the shearing sheet bound- 
ary condition, which allows such nonlinear coherent modes 
to develop. This behavior is not necessarily unphysical or 
irrelevant to actual disk systems, but this point needs to be 
checked in the future. 

Finally, let us come back to the magnetic Prandtl num- 
ber dependence of a. As pointed out earlier, the dependence 
of the transport efficiency on the magnetic Prandtl number 
indicates a back-reaction of small scales on large ones. We 
make here a few comments on this feature. In particular, we 
shall argue that this behavior must saturate at low and large 
enough Pm. The magnetic Prandtl number is related to the 
ratio of the viscous l v and resistive dissipation scales, 
the exact relation depending on the shape of the turbulent 
energy spectrum. Generally speaking, the Prandtl number 
varies monotonically with the ratio l^/l v , and one expects 
Pm <C 1 (resp. Pm 1) when l v /l v <C 1 (resp. S> 1). 

The spectrum of the largest scales tends to be natter than 
usual turbulent spectra due to the role of the linear instabil- 
ity, down to the scale where the magnetic tension prevents 
the instability to occur (most probably, this "instability sec- 
tor" of the spectrum only represents a small part of the 
overall turbulent spectra of actual disks, because of their 
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enormous Reynolds numbers). Leaving aside these largest 
scales, for Pm <C 1, the spectrum is expected to be Kol- 
mogo rovian and anisotropic dow n to the resistive dissipation 
scale dGoldreich fe S ridhar 1995), while below this scale and 
down to the viscous scale, the velocity spectrum is the usual 
Kolmogorov velocity spectrum and the magnetic spectrum 
drops much faster. On the other hand, for Pm 3> 1, the spec- 
trum shoul d be Kolmogorovian down to the viscous dissipa- 
tion scale l|Goldreich fc Sridharlll995T ). while the magnetic 
spectrum should scale like fc -1 below the viscous diss ipa- 
tion scale and down to resistive scale (|Cho et al.ll2003l ). It 
is therefore tempting to see in a difference of accumulation 
of magnetic energy at small scales the cause of the back- 
reaction of these scales to the largest ones, which would 
create the observed magnetic Prandtl number dependence 
of the turbulent transport efficiency. Nevertheless, in both 
small and large Prandtl number settings, turbulent motions 
in the inertial range are random in phase, so that one ex- 
pects that to lowest order, coupling of the turbulent spec- 
trum with the largest MRI unstable scales vanishes. To next 
order, the steepness of the Kolmogorov spectrum indicates 
that the strength of the coupling decreases with increasing 
Reynolds number in the vicinity of the viscous dissipation 
scale, suggesting that at large enough Reynolds numbers, 
the Prandtl dependence should saturate (a potential caveat 
to this argument being the possible role played by a small 
scale field generation through dynamo action) . Such a satu- 
ration was not observed in our simulations, although a weak 
dependence of our results on the magnitude of the Reynolds 
number may be detected on Fig. 1101 however, such an effect 
might also arise from resolution requirements, which makes 
our lower Reynolds number results confined to the larger 
Prandtl number range. Unfortunately, our results can hardly 
be improved upon with the present generation of computers, 
leaving the question of the Reynolds number saturation of 
the Prandtl number dependence open, as well as the over- 
all difference in transport efficiency between the small and 
large Prandtl number cases. Resolving this issue is crucial 
to ascertain the role of the magneto-rotational instability in 
disk transport. 
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